
outcomes <- c("SSDI","SSI","educy","hhinc","emp")


gc()
source("./sipp_eventpreamble.R")


eventdata<-sipp
eventdata[,highedu:=eeducate>39]
eventdata[,preems:=ems[time_limit==0],by=.(spanel,ssuid,epppnum)]
eventdata[preems%in%c("Married, spouse present","Never Married"),premarried := preems == "Married, spouse present"]
eventdata[workcheck==0,premarried := NA]
eventdata[,preobs:=max((time_limit==0)),by=.(spanel,ssuid,epppnum)]

wealthvals <- eventdata[,max(prewealth),by=.(spanel,ssuid,epppnum)]$V1
wealthvals<-quantile(wealthvals[!is.na(wealthvals)])
eventdata[,wealthy:=max(which(max(prewealth) > wealthvals)),by=.(spanel,ssuid,epppnum)]
eventdata[is.infinite(wealthy),wealthy:=NA]
#eventdata[,wealthy:=prewealth > 50000]

eventdata[,preearn:=max(tpearn_real[time_limit<0]),by=.(spanel,ssuid,epppnum)]
preearnvals <- eventdata[ems=="Married, spouse present" & preearn>0,max(preearn,na.rm=TRUE),by=.(spanel,ssuid,epppnum)]$V1
preearnvals<-quantile(preearnvals[!is.na(preearnvals)])
eventdata[,prework:=max(preearn > 0,na.rm=TRUE),by=.(spanel,ssuid,epppnum)]
eventdata[,prewage:=max(which(max(preearn) > preearnvals)),by=.(spanel,ssuid,epppnum)]
eventdata[is.infinite(prewage),prewage:=NA]
eventdata[is.na(premarried),prewage:=NA]
eventdata[,prewage:=prewage+1]
eventdata[premarried==0,prewage:=0]
eventdata[premarried==1 & prework==0,prewage:=1]


sp_earnvals <- eventdata[ems=="Married, spouse present" & sp_tpearn_real>0,max(sp_tpearn_real,na.rm=TRUE),by=.(spanel,ssuid,epppnum)]$V1
sp_earnvals<-quantile(sp_earnvals[!is.na(sp_earnvals)])
eventdata[,sp_work:=max(sp_tpearn_real > 0,na.rm=TRUE),by=.(spanel,ssuid,epppnum)]
eventdata[,spwage:=max(which(max(sp_tpearn_real) > sp_earnvals)),by=.(spanel,ssuid,epppnum)]
eventdata[is.infinite(spwage),spwage:=NA]
eventdata[is.na(ems) | !ems%in%c("Married, spouse present","Never Married"),spwage:=NA]
eventdata[,spwage:=spwage+1]
eventdata[ems%in%c("Never Married"),spwage:=0]
eventdata[ems%in%c("Married, spouse present") & sp_work==0,spwage:=1]


eventdata[,sp_preearn:=max(sp_tpearn_real[time_limit<0]),by=.(spanel,ssuid,epppnum)]
sp_preearnvals <- eventdata[ems=="Married, spouse present" & sp_preearn>0,max(sp_preearn,na.rm=TRUE),by=.(spanel,ssuid,epppnum)]$V1
sp_preearnvals<-quantile(sp_preearnvals[!is.na(sp_preearnvals)])
eventdata[,sp_prework:=max(sp_preearn > 0,na.rm=TRUE),by=.(spanel,ssuid,epppnum)]
eventdata[,sp_prewage:=max(which(max(sp_preearn) > sp_preearnvals)),by=.(spanel,ssuid,epppnum)]
eventdata[is.infinite(sp_prewage),sp_prewage:=NA]
eventdata[is.na(premarried),sp_prewage:=NA]
eventdata[,sp_prewage:=sp_prewage+1]
eventdata[premarried==0,sp_prewage:=0]
eventdata[premarried==1 & sp_prework==0,sp_prewage:=1]

eventdata[,time_limit:=floor(time_limit)]
eventdata[,age_limit:= floor(age_limit)]
eventdata[,time:=as.factor(time_limit)]
eventdata[,time:=relevel(time,ref="-1")]
eventdata[,N:=1]
eventdata[time_limit<=0&age_limit<62 & preobs==0,table(time_limit)]
eventdata[,marriedwealthy:=interaction(married,wealthy)]
eventdata[,spwagewealthy:=interaction(spwage,wealthy)]




write.csv(feols(applyorDI ~ as.factor(time)-1,
                data= eventdata[as.numeric(as.character(time))>=-4 & as.numeric(as.character(time)) <= 4,],
                cluster="ssuid")$coeftable[,c("Estimate","Std. Error")],
          file=paste0("./SIPP_eventstudy/headsev_agebal_level_applyorDI.csv"), 
          row.names = TRUE, col.names = FALSE)


eventdata_orig<-copy(eventdata)

groupvars <- c("married","marhome","spwage","ownhome")
rownums<-c(1,2,2,1)
#"spwagewealthy"
grouprefs <- c(TRUE,"TRUE.TRUE",3,TRUE)
vlabs<-list(c("Single","Partnered"),
            c("Single Renter","Single Owner","Partnered Renter","Partnered Owner"),
            c(#"Single","Non-Working Spouse",
              "1st Quartile","2nd Quartile","3rd Quartile","4th Quartile"),
            c("Renter","Homeowner"))
vlines<-list(c("dashed","solid"),
             c("dashed","dotted","dotdash","solid"),
             1:4,
             c("dashed","solid"))

feregs<-list()
feregs_contrwage<-list()
feregs_uw<-list()
feregs_dynamic<-list()
feregs_trend<-list()
feregs_post1<-list()
for(vnum in 1:length(groupvars)){
  feregs[[vnum]]<-list()
  feregs_contrwage[[vnum]]<-list()
  feregs_uw[[vnum]]<-list()
  feregs_dynamic[[vnum]]<-list()
  feregs_trend[[vnum]]<-list()
  feregs_post1[[vnum]]<-list()
  v<-groupvars[vnum]
  vref <- grouprefs[vnum]
  ##########Balance weights:##############
    eventdata<-eventdata_orig[highedu==0,]
  
  if(v %in%c("spwage")){
    eventdata <- eventdata[get(v) >= 2,]
  }
  vvals <- unique(eventdata[,get(v)])
  vvals<-vvals[!is.na(vvals) & vvals!=vref]
  
  ylabs<-c("Share Receiving DI Benefits",
           "Share Applied for DI Benefits \n in Last 12 Months",
           "Share Received or Applied for DI Benefits \n in Last 12 Months",
           "Share Working Full-Time",
           "Average Annual Earnings")
  outvars<-c("essself","eapplyss","applyorDI","emp_full","tpearn_real")
  onum<-1
  for(out in outvars){
  eventdata[,pweight:=NULL]
  eventdata[,dummy:=get(v)==vref]
  for(vlev in vvals){
    eventdata[,eval(paste0("pval_",vlev)):=NULL]
    eventdata[!is.na(get(v)) & !is.na(get(out)) & !is.na(agebin) & get(v)%in%c(vref,vlev) 
              ,eval(paste0("pval_",vlev)):=feols(as.formula(paste0("dummy~ 1 | interaction(as.factor(time_limit),as.factor(agebin))")),
                    data = eventdata[!is.na(get(v)) & !is.na(get(out)) & !is.na(agebin)  & get(v)%in%c(vref,vlev),
                                     ],)$fitted.values]
  eventdata[get(v)==vlev & get(paste0("pval_",vlev)) < 1 & get(paste0("pval_",vlev)) > 0,
            pweight:=get(paste0("pval_",vlev))/(1-get(paste0("pval_",vlev)))]
  }
  #Imposing common support:
  for(vlev in vvals){
    eventdata[get(v)==vref & !is.na(get(paste0("pval_",vlev))),pweight:=1]
    eventdata[(round(get(paste0("pval_",vlev)),5) == 1 | round(get(paste0("pval_",vlev)),5) == 0),pweight:=NA]
  }
  eventdata[!is.na(pweight),weighted.mean(tage,pweight),by=eval(v)]
  eventdata[,mean(tage),by=eval(v)]
  #######################################
  eventdata[,agebin1:=NULL]
  eventdata[get(v)==1,agebin1:=agebin]
  eventdata[get(v)==0,agebin1:=0]
  
  print(paste0("Unweighted ", v,":", out))
  print(feols(as.formula(paste0(out," ~as.factor(",v,')')),
              data= eventdata[post==1,],
              cluster="ssuid")
  )

  print(paste0("Unweighted by age, ", v,":",out))
  print(summary(feols(as.formula(paste0(out," ~ as.factor(agebin)*as.factor(",v,")- 1")),
              data= eventdata[post==1,],
              cluster="ssuid"))
  )
  
  print(paste0("Weighted ", v,":",out))
  feregs[[vnum]][[onum]] <- feols(as.formula(paste0(out," ~as.factor(",v,')-1')),
                                  data= eventdata[post==1 & as.numeric(as.character(time)) <=4,],
                                  weights=eventdata[post==1 & as.numeric(as.character(time)) <=4,pweight],
                                  cluster="ssuid")

  eventdata[,time0:=as.numeric(as.character(time))==0]
  feregs_post1[[vnum]][[onum]] <- feols(as.formula(paste0(out," ~interaction(as.factor(",v,"),time0) -1")),
                                  data= eventdata[post==1 & as.numeric(as.character(time)) >=0 & as.numeric(as.character(time)) <=4,],
                                  weights=eventdata[post==1 & as.numeric(as.character(time)) >=0 & as.numeric(as.character(time)) <=4,pweight],
                                  cluster="ssuid")
  
  feregs_contrwage[[vnum]][[onum]] <- feols(as.formula(paste0(out," ~as.factor(",v,') + as.factor(prewage) -1')),
                                  data= eventdata[post==1 & as.numeric(as.character(time)) <=4,],
                                  weights=eventdata[post==1 & as.numeric(as.character(time)) <=4,pweight],
                                  cluster="ssuid")
  
  feregs_uw[[vnum]][[onum]] <- feols(as.formula(paste0(out," ~as.factor(",v,')-1')),
                                  data= eventdata[post==1 & as.numeric(as.character(time)) <=4,],
                                  cluster="ssuid")
  print(summary(feols(as.formula(paste0(out," ~as.factor(",v,')')),
                      data= eventdata[post==1,],
                      cluster="ssuid")))
  
  
  meandata<-eventdata[!is.na(get(v)) & !is.na(get(out)),
                      .(outcome=mean(get(out)),N=sum(N)),by=c("time",v)]
  meandata<-meandata[order(time),]
  
  
 pdf(paste0("./SIPP_eventstudy/headsev_level_",out,"_by",groupvars[vnum],".pdf"))
  print(ggplot(data= meandata[as.numeric(as.character(time))>=-4 & as.numeric(as.character(time)) <= 4,])+
          geom_line(aes(x=as.numeric(as.character(time)),y=outcome,
                        group=as.factor(get(v)),
#                        color=as.factor(get(v)),
                        linetype = as.factor(get(v))),size=1)
        +geom_vline(xintercept=0,linetype="dashed")
        +geom_hline(yintercept=0,linetype="dashed")
#       +ylim(c(0,0.5))
#        +scale_x_continuous(breaks=-4:4)
#        +scale_linetype_manual(labels=c("Single","Partnered"),values=c("dashed","solid"))
        +scale_linetype_manual(labels=vlabs[[vnum]],values=vlines[[vnum]])
        +theme(legend.position="bottom")
        +labs(x="Years from Onset",y=ylabs[onum],linetype="",color="")
        +theme(legend.key.size = unit(2,"line"))
  )
  dev.off()
  
  meandata<-eventdata[!is.na(get(v)) & !is.na(get(out)) & !is.na(pweight),
                      .(outcome=weighted.mean(get(out),pweight),N=sum(N)),by=c("time",v)]
  meandata<-meandata[order(time),]

  feresult<-feols(as.formula(paste0(out," ~as.factor(",v,')')),
                  data= eventdata[as.numeric(as.character(time))>=0 & as.numeric(as.character(time)) <= 4,],
                  weights=eventdata[as.numeric(as.character(time))>=0 & as.numeric(as.character(time)) <= 4,pweight],
                  cluster="ssuid")
  summary(feresult)
  
  
  
  feregs_dynamic[[vnum]][[onum]]<-feols(as.formula(paste0(out," ~ as.factor(time) + as.factor(time):",v,'-1')),
              data= eventdata[as.numeric(as.character(time))>=-4 & as.numeric(as.character(time)) <= 4,],
              weights=eventdata[as.numeric(as.character(time))>=-4 & as.numeric(as.character(time)) <= 4,pweight],
              cluster="ssuid")

  feregs_trend[[vnum]][[onum]]<-feols(as.formula(paste0(out," ~ as.numeric(as.character(time)) + as.numeric(as.character(time)):",v,'+ ',v)),
                                        data= eventdata[as.numeric(as.character(time))>=-1 & as.numeric(as.character(time)) <= 4,],
                                        weights=eventdata[as.numeric(as.character(time))>=-1 & as.numeric(as.character(time)) <= 4,pweight],
                                        cluster="ssuid")
  

  pdf(paste0("./SIPP_eventstudy/headsev_agebal_level_",out,"_by",groupvars[vnum],".pdf"))
  print(ggplot(data= meandata[as.numeric(as.character(time))>=-4 & as.numeric(as.character(time)) <=4,])+
          geom_line(aes(x=as.numeric(as.character(time)),y=outcome,
                        group=as.factor(get(v)),
                        linetype = as.factor(get(v))),size=1)
        +geom_vline(xintercept=0,linetype="dashed")
        +geom_hline(yintercept=0,linetype="dashed")
        +scale_linetype_manual(labels=vlabs[[vnum]],values=vlines[[vnum]])
        +theme(legend.position="bottom")
        +labs(x="Years from Onset",y=ylabs[onum],linetype="",color="")
        +theme(legend.key.size = unit(2,"line"))
        +guides(linetype=guide_legend(nrow=rownums[vnum],byrow=TRUE))
  )
  dev.off()
  
  write.csv(feregs_dynamic[[vnum]][[onum]]$coeftable[,c("Estimate","Std. Error")],
            file=paste0("./SIPP_eventstudy/headsev_agebal_level_",out,"_by",groupvars[vnum],".csv"), 
            row.names = TRUE, col.names = FALSE)
  
  
  onum<-onum+1
  }
  names(feregs[[vnum]])<-outvars
  names(feregs_contrwage[[vnum]])<-outvars
  names(feregs_uw[[vnum]])<-outvars
  names(feregs_dynamic[[vnum]])<-outvars
  names(feregs_trend[[vnum]])<-outvars
  names(feregs_post1[[vnum]])<-outvars
}

names(feregs)<-groupvars
names(feregs_contrwage)<-groupvars
names(feregs_uw)<-groupvars
names(feregs_dynamic)<-groupvars
names(feregs_trend)<-groupvars
names(feregs_post1)<-groupvars

feregs_trend$married$essself$coeftable
feregs_trend$ownhome$essself$coeftable

